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ABSTRACT 

To date, no self-consistent numerical simulation of the solar interior has suc- 
ceeded in reproducing the observed thinness of the solar tachocline, and the per- 
sistence of uniform rotation beneath it. Although it is known that the uniform 
rotation can be explained by the presence of a global-scale confined magnetic 
field, numerical simulations have thus far failed to produce any solution where 
such a field remains confined against outward diffusion. We argue that the prob- 
lem lies in the choice of parameters for which these numerical simulations have 
been performed. We construct a simple analytical magneto-hydrodynamic model 
of the solar interior and identify several distinct parameter regimes. For realistic 
solar par ameter values, our results are in broad agreement with the tachocline 



model of iGough fc Mclntyrd . In this regime, meridional flows driven at the base 
of the convection zone are of sufficient amplitude to hold back the interior mag- 
netic field against diffusion. For the parameter values used in existing numerical 
simulations, on the other hand, we find that meridional flows are significantly 
weaker and, we argue, unable to confine the interior field. We propose a method 
for selecting parameter values in future numerical models. 

Subject headings: MHD — Sun: interior — Sun: magnetic fields — Sun: rotation 



1. INTRODUCTION 



The differential rotation of the Sun's convective envelope terminates abruptly at the 
interface with the underlying radiative zone. The transition to uniform rotation takes 
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place across a stably-stratified shear layer known as the solar tach ocline (Spiegel &: Zahn 



19921). whose thickness is at most a few percent of the Sun's radius (IThompson et al. 



Kosovichev et al 



Basu fc Antiall2003[ ) 



19971 ISchou et al.lll998l : Elliott fc Goughlll999l : ICharbonneau et al. 



1996 



19991 : 



More than tw o decades after being first observed (IChristensen-Dalsgaard fc Schoulll988 
Brown et al]ll989h a complete theory of the dynamics of the tachocline is still lacking. In 
particular, the thinness of the tachocline seems to be at odds with the known angular mo- 
mentum tra nsport properties of stratified rota. ting fluids, including advection by meridional 
circulations dSpiegel fc Zdi^ll992l : lElliottlll997l ). Such circulations have a tendency to burrow 
into the radiative zone, transporting the convection zone's angular momentum and thereby 
thickening the tachocline. In order for the tachocline to remain thin, some additional mecha- 
nism must be present that transports angular momentum latitudinally, and in such a way as 
to enforce uniform rotation. In addition, the angular velocity of the radiative zone inferred 
from helioseismology lies within the range of angular velocities observed at th e solar sur- 
face, indicating that the spin-down of the surface by magnetic solar-wind drag (jSchatzman 
19621 ) has been communicated throughout the solar interior. This requires significant vertical 
transport of angular momentum. Again, the transport must be such as to enforce uniform 
rotation. We call such transport "frictional" , meaning down-gradient in angular velocity. 

Several authors have argued that a combination of anisotropic turbulence and internal 
wave bre aking can provide the required latitudinal and vertical angular niomen tum trans- 
port (e.g 



kmg can provide tne requirea iatituamai and vertical angular niomen tum trans- 
Kumar fc Quataert 199?! : Zahn et al. 1997 : Charbonnel &: Talon 2005[). However, 



both of these processes are known to ha ve antz-frictional pr operties (e.g., iMcIntyrd Il994 



Gough fc Mclntyrd Il998l ). On this basis, iGough fc Mclntyrd argued that no non-magnetic 



model can explain the observed thinness of the tachocline and the persistence of uniform ro 
tation within the radiative zone. On the other hand, the presence of a global-scale primordial 



magnetic field provides a natural explana tion for the interior's uniform rotation (e.g.. lFerraro 



19371 : lMestellll953l : iMestel &: Weisall987l ). Such a field enforces uniform rotation through the 
elasticity of the field lines imparted by magnetic tension. The thinness of the tachocline 
can also be explained by the presence of such a field, provided that the field lines are very 



to high latitudes by Maxwell stresses ( 


Riidiger & Kite 


latinov 


1997 


Gough & Mclntvre 


1998; 


MacGregor & Charbonneau 


1999; 


Wood & Mclntvrel 


2011 


) . We call such a field "confined" , 



meaning that the time-averaged mean field resides below the base of the convection zone. 
An unconfined field, on the other hand, would cause the convection zone's differential rota- 



tion to propagate i nto the radiative zone (IMacGregor fc Charbonneaul Il999l : iGaraudI 12002 
Brun fc Zahnll2006h . 
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To explain the confinement of tlie magnetic field, some process must be invoked to 
counteract tlie outward diffusion of tfie field throughout the Sun's lifetime. Global-scale 
meridional flows that downwell from the convection zone into the radiative zone offer a pos- 
sible mechanism for achieving magnetic field confinement — these are the same meridional 
flows that would lead to tachocline thickening in the absence of an interi o r mag netic field. 
This confinement mech anism was first su g gested by iGough fc Mclntyrd (119981 ) and later 
studied numerically by iGaraud fc GaraudI (120081 ). The fiows are generated by gyroscopic 
pumping in the convection zone, a mechanism which has been widely studie d in the context 
of the Earth's atmosphere, and niore recently in an astrophysica l context (iMcIntyrd 120071 : 
Garaud fc Acevedo-ArreguinI l2009l : iGaraud fc Bodenheimerl |2010| ) . In brief, the same angu- 
lar momentum transport by Reynolds stresses that gives rise to the differential rotation of 
the convection zone als o drives a meridional circulation through angular momentum conser- 
vation (see Figure 1 of IGaraud fc Bodenheimerl |2010| ). Such meridional circulations are a 
robust consequence of anisotropic angular momentum transport, as originally discussed by 
KippenhahnI (119631 ). 



Although the Sun's meridional flows can be observed directly at the solar surface 



and inferred in the near- s urface layers frorn local helioseismology (e.g., iHaber et al.l 12002 



Zhao &: KosovichevI |2004J : iGizon et al.l |2010| ) , their amplitude and structure deeper within 
the solar interior is currently unknown. Theoretical and numerical models of the large-scale 
dynamics of the solar interior are required in order to establish whether the gyroscopically 
pumped circulation is strong enough, and penetrates deeply enough into the radiative zone, 
to conflne the interior magnetic fleld and thereby explain th e observed ta c hocline struc- 
ture. Several such nurnerical studies have been p erformed ( iGaraudI l2002t iBrun fc Zahn 
20061 : IGaraud fc GaraudI l2008l : IStrugarek et al.ll2011 ). but none has successfully reproduced 
the fleld- conflnement scenario of iGough fc Mclntyrd . However, since computing limitations 
prevent any numerical model from reaching true solar parameter values, this failure may sim- 
ply reflect the fact that the numerical models are not in the parameter regime of relevance 
to the solar tachocline. 



An alternative, analytical approach was proposed by IGaraud &: Brummelll (|2008[ ). and 
further developed by IGaraud &: Acevedo-ArreguinI (120091 . hereafter IGAA09I ). These prelim- 
inary models examined the amplitude and depth of penetration of global-scale meridional 
flows into the radiative zone in the absence of an interior magnetic fleld. In these non- 
magnetic models, the burrowing of the meridional flows was halted only by the presence of 
viscosity. In the steady state, the meridional mass flux was found to decay exponentially 
with depth below the radiative-convective interface on a lengthscale ~ Rq/o', where Rq is 
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the solar radius and 



N 



Here u is the viscosity, k, is the thermal diffusivity, N is the buoyancy frequency, and Qq is the 
mean solar rotation rate. This result shows how the burrowing tendency of the meridional 
flows is strengthened by rotation, and weakened by stable stratification and viscosity. Within 
the solar tachocline a < 1, and so viscosity alone cannot prevent the meridional flows from 
burrowing into the radiative zone. 



Garaud &: Brummelll (|2008[ ) found that the amplitude of the steady-state meridional 
flows in their model was sensitive not only to the value of a, but also to conditions at th e 
interface between the c onvective and radiative zones ( see also iBretherton fc Spiegellll968l ). 
In a subsequent study, iGaraud fc Bodenheimerl ( I2OIOI ) showed that the vertical mass flux 
below the radiative-convective interface can be quantified in terms of two constraints. The 
first constraint is that there must be stresses present in the radiative zone that overcome 
Taylor-Proudman balance. In the absence of such stresses, any meridional flows would 
have to be parallel to the rotation axis, which is not compatible with mass conservation. 
When such stresses are present, however, downwelling meridional flows are able to turn 
around and return to the convection zone — this is another exarnple of gyroscopic pumping. 
Garaud &: Bodenheimerl called this the "mechanical" constraint]^ 



The second constraint is due to the presenc e of stable stratification, whic h inhibits 
vertica l flows; we call this the "thermal" constraint. iGaraud fc Bodenheimerl ( I2OIOI . hereafter 
GBlOl ) found that this places an upper bound on the vertical flow velocity W of the form 



W < 



Gt 



ES 



where tss is the local Eddington-Sweet timescale. 



(2) 



(3) 



[e.g., iSpiegel fc Zahnlll992l ). and G is a dimensionless "geometrical" factor that depends on 
the details of the star's internal structure. 

The two constraints set two upper bounds for the meridional flow velocity, which then 
scales as the smaller of the two. In the case where the mechanical constraint is smaller, 



^ Since the mechanical constraint arises from the need to balanc e azimuthal forc es, it applies only to the 
steady state. Transient meridional flows, such as those studied by ISpiegel &: ZahnL are not subject to this 
constraint and can be much stronger. 
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the meridional flow veloci ty depends on the nature of the stresses that overcome Taylor- 



Proudman balance. IGBIOI considered both laminar viscous stresses and turbulent Reynolds 
stresses in the convective cores of "lithium-dip" stars. In the present paper we consider 
magnetic stresses, which we argue are the dominant mechanism for angular momentum 
transport in the solar radiative zone. For simplicity, we neglect any turbulent stresses outside 
of the convection zone. 



Following iGAAOQl and IGBIOI we use a linearized, steady-state, Cartesian model of the 
solar interior to study the interaction between meridional flows and an imposed conflned 
magnetic fleld. In ||2] we describe the Cartesian model and derive the governing equations. 
In ^we present numerical and analytical solutions in the case where stratiflcation is absent. 
We demonstrate a simple relationship between the structure and magnitude of the confined 
magnetic fleld and the pattern and amplitude of the meridional flows. In particular, we 
quantify the mechanical constraint from the conflned magnetic fleld. In ^we extend these 
results to the stratifled case, recovering the thermal constraint mentioned above. We then 
discuss in ^ what insight the results of our simplifled model provide into the nonlinear 
dynamics of the solar interior. In particular, we identify the various parameter regimes that 
occur in our model, and we discuss in which of these regimes the predicted meridional flows 
could plausibly bring about conflnement of the magnetic fleld. The results show that great 
caution must be exerted when interpreting the results of numerical simulations that use non- 
solar values for the governing parameters, and provide guidance as to the design of future 
numerical simulations of the solar interior. 



A Cartesian model 



We use a similar model to that of IGAAOQI . Although our model is intended as a sirt i- 
plifled Cartesian analogue of the tachocline picture proposed by lGough fc Mclntyrd (Il998l ). 
it also goes beyond theirs in certain respects. We explicitly model the thermal and dynam- 
ical coupling between the convective and radiative zones, and are thus able to draw more 
quantitative conclusions concerning the role of this interaction in the global-scale solution. 
We also consider much wider numerical ranges in the governing parameters, and describe 
the structure of the solutions in several distinct parameter regimes. Our aim is to identify 
in which parameter regime, if any, the propagation of the differential rotation below the 
convection zone is restricted to a thin layer with a tachocline-like structure. 
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2.1. Model geometry 

We measure distances in units of the solar radius, Rq ~ 7 x lO^^cni, and velocities in 
units of RqQq, where Qq ~ 3 x 10~^s~^ is the mean solar rotation rate, which we use as the 
inverse timescale. As illustrated in figuredl the Cartesian coordinate system [x, y, z) is chosen 
such that X is the azimuthal coordinate, y G [0, tt] is the latitudinal coordinate, and z G [0, 1] 
is the vertical coordinate. Although there is no precise equivalence between our Cartesian 
geometry and the spherical geometry of the solar interior, for the sake of interpretation we 
will take ?/ = and y = vr to be the "poles", and ?/ = 7r/2 to be the "equator". The 
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Fig. 1. — A meridional cross-section through the Cartesian modeL The y and z axes represent latitude 
and altitude respectively. The shaded area marks the convection zone, where a forcing is applied to mimic 
the generation of differential rotation by turbulent stresses. The rotation axis is vertical. A global-scale 
magnetic field Bo(z) ~ Boe^^^^ey is confined by uniform downwelling Uq. 



radiative-convective interface is at z = /i ~ 0.7, so that h < z < 1 represents the convection 
zone, and < z < h represents the radiative zone. 



2.2. Background state 

The steady background state considered is depicted in figure [H We model the transition 
between the convective and radiative zones by imposing a vertical profile for the dimensionless 
buoyancy frequency n{z) with 

n2(z) = ^(l + tanhf^')l , (4) 



2 V A 
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so that 77, = in the convection zone and n = rirz > in the radiative zone. The lengthscale 
A is introduced in order to ensure smoothness of the background state, which is necessary 
for numerical reasons. We may regard A as the thickness of the convective overshoot region 
at the base of the convection zone. 

We impose a horizontal background magnetic field Bq of the form 

Bo = Bo{z)ey, (5) 

where ey is the unit vector in the y direction. In order to confine this field against ohmic 
diffusion, in the steady state, we must invoke a background meridional flow. The simplest 
case is that of uniform downwelling of dimensionless magnitude Uq, and 

Bo{z) = B,e-''\ (6) 

where Bq is a constant and 5 is the dimensionless lengthscale 

Here rj and represent the magnetic diffusivity in dimensional and nondimensional units 
respectively; we call the "magnetic Ekman number". Although this background state 
cannot hold over all latitudes and depths within the solar interior (for reasons of mass conser- 
vation and solenoidality of Bq), we hope that it offers a reasonable qualitative approximation 
to conditions within the tachocline at middle and high latitudes. 

We note that Bq is not a force-free field since it has a non-constant magnetic pres- 
sure i?o/(87r). However, the magnetic pressure can be balanced by a perturbation to the 
background gas pressure po- Estimates of the str ength of the Sun's interior magnetic field 



typic ally lie within the range 10 '^-10 gauss (e.g. iMestel fc Weiss! 119871 : iGough fc Mclntyre 



19981 . see also §4.3p . For field strengths in this range, the gas pressure far exceeds the 



magnetic pressure, and so the required perturbation is very small. 

In what follows we consider linear, steady-state perturbations to this background state. 
The principle shortcoming of our model is, arguably, the artificial manner in which the mag- 
netic field is confined to the interior, via the imposed downwelling Uq. Magnetic confinement 
is thereby built into the model, and controlled by the parameters Bq and 6. This idealiza- 
tion is necessary to allow a linear study. A complete model of the tachocline would need 
to describe self-consistently the processes that act to confine the magnetic field, rather than 
treating it as part of the background. Such a model would necessarily be nonlinear, and 
awaits future numerical work. We note also that our Cartesian model cannot adequately 
describe the polar regions, where the effects of geometrical curvature beco me significant. The 



polar magnetic confinement problem has been studied in detail elsewhere ( jWood fc Mclntyre 



20111 ). and self-consistent, fully nonlinear solutions have been obtained. 
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2.3. Model equations 



As in tlie studies of IGAA09I and lGBlOl we do not model the turbulence in the convection 
zone in detail. Instead, we introduce a simple parametric model to describe its role in driving 
a large-scale differential rotation profile, and in enhancing the transport of heat. In the 
region z G [/i, 1] we model the turbulent heat transport as a diffusive process, and we replace 
the divergence of the Reynolds stress with a forcing term that represents linear relaxation 
towards a prescribed rotation profile, which we choose to mimic the observed differential 
rotation of the solar convection zone. 



For simplicity we also make the Boussinesq approximation (e.g. lSpiegel fc Veronislll960l ). 
in which we neglect variations in the background pressure po? density po and temperature 
To, and we neglect pressure perturbations in the equation of state. Although the Boussinesq 
approximation cannot be rigorously justified over the entirety of our domain, we do not 
expect this approximation to have a qualitative effect on our results. Indeed, we are interested 
primarily in solutions for which the meridional flows and differential rotation are conflned 
to a thin, tachocline-like layer below the convection zone, within which the Boussinesq 
approximation should be quite accurate. 

The nondimensional governing equations, linearized about the background state de- 
scribed in the previous section, are 

'1 



2e^ X u = -Vp + Te^ + Ae"^/"^ 



■e^: X B + (V X B) X e„ 



n^(z)w = E,V-(/(^)VT) (9) 

= V X (u X e^^/^e^ - ^e^ x B - V x B) (10) 

V-u = (11) 
V-B = 0, (12) 

where u = [ux-, Uy, Uz) is the velocity perturbation, B = {hx, by, hz) is the perturbation to the 
magnetic fleld, p is the pressure perturbation, and T is the temperature perturbation. These 
equations have been scaled using pqRqQq as the unit of pressure, TqRqVI?^/ g as the unit of 
temperature, and i3o/E^ as the unit of magnetic fleld strength. The following dimensionless 
parameters have been introduced: 

A = — , the Elsasser number at ^ = 0; (13) 

Ey = , the Ekman number; (14) 

~ d2 o ' thermal Ekman number in the radiative interior. (15) 
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Note that the non-dimensionahzed perturbation equations have no exphcit dependence on 
the magnetic Ekman number E^. However, because has been used to scale the magnetic 
field, its value governs the magnitude of perturbations relative to the background field Bq. 

The Elsasser number A is a measure of the relative magnitude of Lorentz forces compared 
to Coriolis forces. It is convenient to define also a local Elsasser number, 

Aioc(.) = (16) 

= Ae-2^/^ (17) 
in order to measure the relative magnitudes of these forces at each altitude z. 

In both the momentum equation ([8]) and the thermal energy equation (Q we have ne- 
glected terms describing advection by the background downwelling Uq. It can be shown that, 
in the results presented here, including these terms would produce only a small correction 
to the solutions. Neglecting these terms simplifies the analysis and reduces the dependence 
of the solutions on the artificially imposed fiow Uq. 

The final term in the linearized momentum equation ([8]) describes a linear relaxation 
towards a prescribed rotation profile Ucz. The dimensionless parameter A, which determines 
the rate of relaxation, can be regarded as the inverse of the convective turnover timescale. 
This parame terization for the angu l ar mo r nentum t ran sport b y Reynolds stresses is similar to 



that used by lBretherton &: Spiegell (ll968[ ). lGAA09l and lGBlOl . For simplicity, we assume that 



A takes the constant value 1/tc in the convecti on zone, and that the prescribed differential 



rotation u^z is sinusoidal in latitude. Following IGAA09I . we adopt the profiles 

Ucz = Mcz(2;)e'''^e^, (18) 

and A(;z) = ^|l + tanh(^^^^| , (19) 

where A is the same "overshoot" length as in equation (jlj). In the numerical results presented 
later we set Tc = 0.1 and k = 2, representing equatorial symmetry at the solar surface, and 
Ucz{z) is taken to be linear in z. We show in appendix |B] that the fiow below the convection 
zone is not very sensitive to the value of Tc, or to the form of Ucz{z), and that the results are 
easily extended to any particular Ucz{z). This insensitivity to the choice of parameterization 
within the convection zone refiects the robustness of both the gyroscopic pumping mechanism 
and the burrowing tendency of meridional fiows. 

Finally, we model turbulent heat transport within the convection zone as an increase in 
the thermal diffusivity by a factor / relative to its laminar, microscopic value. We take 

m = 1 + ^ {l + tanh (^) I , (20) 



so that / = 1 within the radiative zone and / = /o ^ 1 within the convection zone. Ther- 
mal perturbations at the radiative-convective interface are therefore communicated more 
efficiently into the convection zone than into the radiative zone. In the limit /o — )■ oo we 
expect the convection zone to become isothermal, which is equivalent to isentropic in our 
Boussinesq model. In the numerical results presented in ^we take /o = 10^. The physical 
significance of /o, and its effect on the solutions, is described in §4.3[ below equation (!72|) . 

We assume that the viscosity v and the magnetic diffusivity ?7 retain their laminar, 
microscopic values throughout the solar interior. Introducing larger "turbulent" values for v 
and 1] within the convection zone would not significantly affect the results we present here. 

The perturbation fields must have the same sinusoidal dependence on y as the forcing 
term A(z) u^z in (|H]), and so we seek solutions of the form 



The perturbation equations fE])-f lT^ then become a system of ordinary differential equations 
for the perturbation amplitudes q: 




(21) 



- 2m, 




'X 



(22) 




-ikp - -i.e'^l^ + E, 




- k'^u, - \{z)u^ (24) 



(23) 



k%^ 



'x 



(25) 



k%. 



(26) 



(27) 



(28) 



(29) 



The system of equations fl22|) - (!29|) is solved numerically using a two-point boundary-value 
solver, given the boundary conditions listed below. 
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2.4. Boundary conditions 



The upper and lower boundaries are modeled as impenetrable, (viscous) stress-free, 
isothermal and electrically insulating. That is, we impose 



Uz 

ditx 



^ 





dz 
T 






X = 
db 



> at z = and z = 1, 



(30) 



— — = kbz at z = 0, 
dz 



and 



dbz^ 
dz 



—kb, at z = 1. 



(31) 
(32) 



Most of these boundary conditions have been selected for simplicity and convenience. In 
fact, the solutions turn out to be rather insensitive to the choice of boundary conditions in 
the parameter regime of interest, as we have verified by comparing solutions obtained with 
different sets of boundary conditions. 



3. The Unstratified Case 

3.1. Numerical exploration 

In order to delineate the various physical effects present in the solutions, we begin by 
considering the unstratified case, = 0. Effects associated with stratification will be 
described in detail in §11 

A series of results is shown in figure [2l In the absence of any r nagnetic field (i.e.. 



setting A = 0) we recover the non-magnetic, unstratified solution of iGAAOQl (see their 
figure 2). Since the lower boundary is stress-free, the differential rotation imposed within 
the convection zone is able to spread throughout the radiative zone, and establishes a state 
close to Taylor-Proudman balance. In this steady state the only meridional flows are those 
maintained by viscosity, and their magnitude scales with the Ekman number, E,^. The 
meridional flows therefore become vanishingly weak in the limit Ej, — )■ 0. When A > 0, on 
the other hand, the Lorentz force produces departures from Taylor-Proudman balance, and 
maintains a meridional flow even in the limit E^, — > (figure [2]d). 
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Fig. 2. — Profiles of vertical velocity Uz for a series of unstratified solutions (riiz = 0) with different 
magnetic field strengths and Ekman numbers. All the solutions have Ucz{z) = 1 + {z ~ h) in p^ . and 
A = 0.001 and Tc = 0.1 in (|19l) . The magnetic scaleheight, as defined by (O, is (5 = 0.02. Each panel shows 
the solutions for = 10^^ (solid line), E,y = 10^^ (dashed line) and Ejy = 10^'* (dotted line). 
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For A of order unity, the Lorentz forces are most significant close to the lower boundary 
z = 0, within a boundary layer analogous to an Ekman-Hartmann layer. The structure of the 
boundary layer depends on the particular choice of boundary conditions at z = 0, including 
the magnetic boundary conditions. However, we have verified that the meridional flow 
outside the boundary layer is hardly affected by the choice of magnetic boundary conditions. 

For A ^ 1 the radiative zone divides into two dynamically distinct regions, separated 
by a thin internal boundary layer at z = zq, say (see figure The region above the 

boundary layer is close to Taylor-Proudman balance, with all components of u independent 
of z, indicating that Lorentz forces are locally negligible (see figure |8^ in §4.1l for plots of Ux)- 
In the region below the boundary layer, on the other hand, the Lorentz force is dominant, 
and all components of u decay exponentially with depth. From here on we refer to these 
two regions as being "weakly magnetic" and "magnetically dominated" respectively. The 
boundary layer separating these two regions we call the "magnetic transition layer". Since 
the flow within the magnetically dominated region decays exponentially with depth, the 
solution is insensitive to the choice of boundary conditions at 2; = 0. 

The structure of the weakly magnetic and magnetically dominated regions is illustrated 
further in figure [3l which shows the meridional streamlines and magnetic field lines, as well 
as contours of the azimuthal components of the velocity and magnetic fields, in a verti- 
cal cross-section through a typical solution. Within the convection zone, the forcing that 
drives differential rotation also gyroscopically pumps a meridional flow. Below the base of 
the convection zone the forcing is switched off, and Taylor-Proudman balance is achieved. 
As a result, the meridional streamlines and angular velocity contours become aligned with 
the rotation axis, and a portion of the convection zone's differential rotation extends into 
the radiative zone. This differential rotation winds up the magnetic field lines, producing a 
Lorentz torque whose amplitude increases exponentially with depth. Within the magnetic 
transition layer (here around zq = 0.2) the Lorentz torque is strong enough to overcome 
Taylor-Proudman balance, and gyroscopically pumps a meridional flow. The magnetic field 
lines are dragged upward where the flow is upwelling, < y < 7r/2, and pushed down- 
ward where the flow is downwelling, n /2 < y < tt (not shown in the figure). Below the 
transition layer, meridional flows and differential rotation are both strongly suppressed, and 
consequently perturbations to the magnetic field lines are smaller. 

Figure|l]shows the result of increasing A still further. Over many orders of magnitude the 
only effect of increasing A is to raise the position of the magnetic transition layer, shrinking 
the vertical extent of the weakly magnetic region. The dynamics within the convection 
zone remain unaffected until the point at which the transition layer meets the base of the 
convection zone. 
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y y 

Fig. 3. — Meridional cross-section through the solution in figure [2] with A = 10^^ and Ei, = 10~*. The 
left panel shows the azimuthal velocity in color, overlaid with selected meridional streamlines (solid lines 
for anticlockwise flow, dotted lines for clockwise flow). The right panel shows lines of the poloidal magnetic 
field, including the background field Bq, assuming a magnetic Ekman number of = 10"'*. Although field 
lines have been drawn at regular vertical intervals, the magnetic field strength decays exponentially with 
height. The color scale shows bx{x,y)/BQ. 

It is obviously of interest to identify the factors determining the location of the magnetic 
transition layer. We anticipate that this transition occurs where the local Elsasser number 
Aloe, defined in f|T7|l . takes a critical value. We therefore expect the vertical position of the 
transition layer, Zq, to follow a law of the form 

Zq = const. + i(51nA (33) 

as A is increased, where the constant depends on the aforementioned critical value of Aioc. 
It is readily verified from figure H] that zq does indeed follow such a law. In the following 
sections we analyze the structure of the solutions in detail in order to determine both the 
critical value of Aioc and the amplitude of the gyroscopically pumped meridional flow below 
the convection zone. 



3.2. Analytical solution 



In order to understand the results presented in §3.11 more quantitatively, we now derive 
approximate analytical solutions of the governing equations. In §3.2.1l we derive a boundary- 
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layer solution for the magnetic transition layer. In §3.2.2l we use that boundary- layer solution 
to construct a global analytical solution. 



3.2.1. Transition layer solution 

We seek to describe the transition layer between the weakly magnetic and magnetically 
dominated regions of the flow. Since this transition occurs below the base of the convection 
zone we may neglect the terms in the momentum equation involving A. Based on the 
numerical solutions shown above, we anticipate that the thickness of the transition layer is 
of the same order as 5, the scaleheight of the background magnetic field. Since 5 ^ 1/A;, we 
make the boundary-layer approximation 

d^ A d^ 

dz^ / dz^ ^ ^ 

We also neglect viscosity, i.e., we set E^, = 0. It can be verified a posteriori that viscosity 
plays no significant role in the dynamics of the transition layer provided that 5 is much larger 
than the Ekman length, that is, provided that 

5 > ¥}J\ (35) 

This condition is easily satisfied in the solar interior, where Ej^ ~ 10~^^. 

With these approximations, the governing equations f l2^ -f l2^ reduce to the following 

set: 



— 2Uy 




(36) 




= -ikp-H,e-'/^ 




(37) 


dp 
di 


d dz 


(38) 





dz dz"^ 


(39) 





dz dz"^ 


(40) 


■kUy + — — 

dz 


= 


(41) 


ikby + — — 
dz 


= 0. 


(42) 
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We show in appendix |X] that these can be combined into a single equation for -u^: 



d _ 2 

d^ 6 



e 



u, = 0. (43) 



From inspection of (1431) we deduce that the vertical position of the transition layer, zq, is 
given by 

(e^o/V^)^ = iA;4A2 

^Zo = 6\n [k5/V2^ + ^5 In A , (44) 

which is consistent with (1331) . Furthermore, this expression shows that the transition layer 
is located where the local Elsasser number Aioc, defined by ( JT7I) . takes the critical value 
2/(^5)2. This corresponds to a critical magnetic field strength 

5„,..»Ii.40(AV'7^^G. (45) 



k5 \V& J \ ^ 

where po — 0.2gcm~^ and rjQ ^ 400cm2s~^. 

We seek the solution to (1431) that has the correct behavior in both the weakly magnetic 
region, for z — zq 5, and the magnetically dominated region, for z — zq <^ —5. In 
particular, the solution must match onto the differential rotation, = Ut say, in the region 
immediately above the transition layer. Here, and subsequently, we use a subscript "t" to 
refer to the region immediately above the transition layer, which is also the bottom of the 
weakly magnetic region. We show in appendix |X] that the unique solution with the required 
properties is 

= ikdut [/i(C) - f Re jexp (-^C) }] , (46) 
where C = exp((zo — z)/5), and Ji is the integral 



oo 



s ds 



The weakly magnetic and magnetically dominated regions correspond to C ^ 1 and C ^ 1 
respectively. 

Figure shows the vertical profiles of Uz and Ux-, in the vicinity of the magnetic transition 
layer, from the same numerical solution shown in figure [31 Also shown are their analytical 
counterparts, fH6|) and the corresponding analytical profile of Ux- 



Within the magnetically dominated region, where C ^ 1, the flow is exponentially weak. 
From (jinD we find that 

Uz ^ = ikdut exp I I . (48) 
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Within the weakly magnetic region, where ( ^, Uz tends to a constant, Wt say, as expected 
from the Taylor-Proudman constraint. Specifically, 

Wt = Uz\t^=o = -if^*^ ^t- (49) 

The vertical mass flux within the weakly magnetic region is therefore tied to the differential 
rotation via the magnetic transition layer. 



3.2.2. Global solution 

We now use the relation fH9|) derived from our transition-layer solution to construct 
an approximate analytical solution for the global flow. In this way we derive an explicit 
relation between the constant vertical velocity Wt within the weakly magnetic region and the 
prescribed differential rotation Ucz{z) within the convection zone. To construct the global 
solution, we must also match the flow within the weakly magnetic region, zo < z < to 
the flow wit hin the convection zone, h < z < 1. The matching procedure is similar to that 



employed by iGBlOl . and is described in appendix |Bl Here we present only the result for Wx,. 



In the absence of stratiflcation, rirz = 0, we flnd that Wt is given by 



coth 



\ d } TikS 

where d is the "convective" lengthscal^ 



d^^^^ (51) 



k 

and where Ucz is a weighted average of the forcing in the convection zone 

dz Uc7,{z) cosh 
d^; cosh 



Ucz = ""^ „i — d^ ^ j^g2) 

' ' \ — z^ 



h 



d 



Since d > 1/k > {1 — h), the weight distribution in ( l52l) is roughly uniform, and so Ucz is 
numerically close to the vertical average of Ucz{z) over the convection zone. 



^The lengthscale d in this paper is equivalent to the lengthscale ^out in 



GBIO. 
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Figure [6] shows the vertical profiles of Uz in a series of numerical solutions with various 
magnetic scaleheights 6. The bottom panel in figure [6] compares the value of Uz in the weakly 
magnetic region to the value for Wt predicted by flHIl) . The fit is excellent provided that the 
magnetic transition layer is located sufficiently far above the lower boundary, z = 0, and 
sufficiently far below the base of the convection zone, z = h. Using (1441) . these two conditions 
can be expressed as lower and upper bounds on A: 

2 2e2V5 ^ ^ 



3.3. Summary and discussion for the unstratified case 

We have found that a confined magnetic field is able to halt the burrowing of the 
convection zone's differential rotation into the radiative interior, provided that the magnetic 
field strength exceeds a critical value -Bcrit given by fH5|) . Because the field strength in our 
model increases monotonically with depth, the radiative zone divides into an upper, weakly 
magnetic region, wherein Bo{z) ^ -Bcrit, and a lower, magnetically dominated region, wherein 
Bq{z) ^ -Bcrit- In the weakly magnetic region the flow is subject to the Taylor-Proudman 
constraint, and the differential rotation is independent of z. In the magnetically dominated 
region the flow is subject to the Ferraro constraint, which suppresses both differential rotation 
and meridional flows. In our model Bq{z) increases exponentially with depth on a fixed 
lengthscale S, and so the transition between the weakly magnetic and magnetically dominated 
regions occurs across a layer of thickness ~ 6, within which Bq -Bcrit- 

The transition layer regulates the mass flux that downwells from the convection zone, as 
expressed by equation (H9l) . This is closely related to the mechanical constraint mentioned in 
^ and can be explained as follows. The latitudinal differential rotation below the base of the 
convection zone extends downward until it meets the transition layer, where it winds up the 
magnetic field lines, producing a Lorentz torque that overcomes Taylor-Proudman balance 
and gyroscopically pumps a meridional flow. In the steady state, the downwelling mass flux 
within the weakly magnetic region must be exactly that demanded by the transition layer, 
which is expressed by ( l49l) . An analogy can be made with the more familiar Ekman pumping 
that occurs in Ekman layers. Indeed, if the transition layer were replaced by an artificial 
frictional, impenetrable and non-magnetic horizontal boundary, then (H9l) would be replaced 
by the Ekman-pumping formula, 

Wt = --ik6EUt (54) 



(IGAA09I ). where Se = is the nondimensional Ekman- layer thickness. In the solutions 
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presented here we have S ^ 6e, and so the magnetic transition layer pumps much stronger 
meridional flows than would be pumped by an artificial Ekman layer. 

The sol utions we have found bea r many similarities to the tachocline model originally 
proposed by lGough fc Mclntyrd (119981 ). In particular, we have a differentially rotating region 
below the base of the convection zone, which we might call the tachoc line, and a thin magnetic 
boundary layer at the base of this region, which iGough &: Mclntyrd called the "tachopause" . 
The entire region below the tachopause is held in uniform rotation by the confined magnetic 
field. However, an obvious shortcoming of our unstratified soluti ons is that the "tach ocline" 
has no vertical shear, quite unlike the solar tachocline. In the iGough fc Mclntyrd model, 
the tachocline's vertical shear is explained by the presence of stable stratification, via the 
thermal- wind relation. In the following sections we reintroduce stratification into our model, 
in order to better approximate conditions within the solar tachocline. 



4. The Stratified Case 

The importance of stratification in our model can be measured in terms of the dimen- 
sionless parameter n^^/E^, which appears on the left-hand side of ( l27|l . In fact, it follows 
from equation ([2]) that 

^ = ^etES , (55) 

so n^^/E^ is precisely the dimensionless local Eddington-Sweet timescale. In what follows 
we also use a = nj.z(E,^/'E,^Y^'^ as a measure of the stratification. As mentioned in this 
parameter is particularly relevant in cases where buoyancy and viscosity both contribute to 
the balance of forces. All the numerical solutions presented in this section have = 10~^, 
and so the parameters n^^/^K and a contain equivalent information. 



4.1. Numerical exploration 

We begin again by exploring the behavior of the solutions of fl22l)- (!29l) with varying 
governing parameters. In figure [7| we show the vertical profiles of Uz in a series of numerical 
solutions with increasing stratification, measured in terms of a. When a is sufficiently 
small, < 10~^, the profiles are indistinguishable from one another and from the unstratified 
solution, which has a = 0. The same is true for the profiles of Ux in these solutions, which 
are shown in figure [8^. In particular, the weakly magnetic region {zq < z < h) remains in 
Taylor-Proudman balance. We refer to all such cases as "unstratified" , since any effects due 
to stratification are negligible. 
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As cr is increased beyond 10~^ we observe a reduction in Uz in the weakly magnetic 
region (see figure [7]). At the same time, the jump in across the magnetic transition layer 
drops almost to zero, and instead Ux{z) rises smoothly and monotonically between zq and h, 
indicating that thermal-wind balance has been established within the weakly magnetic re- 



gion. For brevity, and because of the obvious similarity with the model of lGough fc Mclntyre 



(Il998l ). from here on we refer to the weakly magnetic region and magnetic transition layer 



in such cases as the "tachocline" and "tachopause" respectively. 

We find that cases with significant stratification can be further categorized as being 
either "weakly stratified", "moderately stratified", or "strongly stratified". These three 
regimes are now described in turn. 



4- 1.1. The strongly stratified regime 



In common with iGaraud fc Brummelll (120081 ) and IGAA09I . we find that the strength of 
the meridional flow Uz decays exponentially with depth below the convection zone on the 
vertical scale l/{ka). If a is sufficiently large, > 1, then this vertical scale is smaller than 
the tachocline thickness, and so meridional flows gyroscopically pumped at the base of the 
convection zone are halted by a combination of buoyancy a nd viscou s forces before reaching 
the tachopause. Following the terminology introduced by IGAAOQI we refer to such cases 
as "strongly stratified". Even with strong stratification, the differential rotation driven 
in the convection zone still spreads into the radiative interior (figure [8^), but by viscous 
diffusion rather than by the burrowing of meridional circulations. The differential rotation 
winds up the magnetic field lines much as in the unstratified cases described in §3], producing 
a Lorentz torque within the tachopause that gyroscopically pumps a meridional flow. Part of 
that meridional flow is launched upward into the tachocline, where its strength decays with 
height on the lengthscale l/{ka). The strongly stratified solutions therefore have meridional 
circulations within two distinct horizontal layers, at the top and bottom of the tachocline. 



4.1.2. The weakly stratified and moderately stratified regimes 

Figure [8]d shows vertical profiles of T from the same series of solutions shown in figure [71 
In all cases with a < 10~^ we find that there is little variation in the temperature field T 
or its vertical gradient dT/dz across the tachopause. However, in cases with a > 10"^ we 
observe a significant change in the temperature gradient ai z = Zq, indicating that the flow 
within the tachopause contributes to the global-scale thermal equilibrium. We refer to cases 
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with 10 ^ ^ cr ^ 10 ^ as "weakly stratified", and cases with 10 ^ ^ o" ^ 1 as "moderately 
stratified" . 

Figure M illustrates how the mass fiux in the tachocline varies with the strength of the 
interior magnetic field in the weakly and moderately stratified regimes. As in the unstrat- 
ified case, we find that the position of the magnetic transition layer moves upwards as A 
increases, in a manner which is more-or-less independent of a. However, in contrast with 
the unstratified case, the strength of the meridional fiow within the tachocline varies with 
A (cf. figure H]). For increasing values of A, the meridional fiow becomes stronger as the 
tachocline becomes thinner. This result has implications for the solar tachocline, as will be 
discussed in §4.2.21 

Figure [TU] shows meridional cross-sections through two solutions in the weakly and 
strongly stratified regimes. As mentioned in ^ the solar tachocline has a ~ 0.2-0.4 



(jGAAOQl ). and therefore is not strongly stratified in the sense used here. We anticipate 
that either the weakly stratified or the moderately stratified regime offers the best approx- 
imation to the dynamics of the tachocline. In the following sections we analyze these two 
regimes in detail in order to quantify more precisely the parameter ranges over which they 
apply. 



4.2. Analytical solution 

We proceed, as in §3.21 by deriving an analytical boundary-layer solution for the tachopause, 
which can then be used to construct a global analytical solution. In order to simplify the 
analysis we will neglect the viscous terms in the equations, and so our analytical solutions 
cannot be applied in the strongly stratified regime. 



4-2.1. The tachopause 

In the unstratified, weakly stratified and moderately stratified solutions presented in 
figure[H)D, the temperature field T{z) is roughly constant in the region z — zq = 0{6). (Recall 
that the moderately stratified cases are characterized by a jump in dT/dz, rather than T.) 
We can therefore analyze the structure of the tachopause in all three regimes simultaneously 
by approximating T as a constant, Tt. As in §3.2.11 we neglect viscosity and make the 
boundary-layer approximation fl51|) . We then recover equations fl5^ - fH21) exactly, except 
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that fl38|) is replaced by 



dp 



5 dz 



(56) 



Within the tachocline, where z — zq ^ 6, all the terms involving A are negligible, and so we 
have 



2u^ 



d~z~^' 



ikuy + 



duz 
dz 





—ikp 


0. 



(57) 
(58) 

(59) 
(60) 



We deduce that 



Ux ~ Ut- ^ikTt{z - zo 

Uy ^ 
Uz Wt 



as 



6 



+00 



(61) 



where the constants Ut and Wt represent the values of itx and Uz immediately above the 
tachopause, or equivalently, at the bottom of the tachocline. 



As in §3.2.11 we can combine all of our boundary-layer equations into a single equation 
for u,: 



f d 2 

Vdi ~ 5 



dz 



+ ifc^A^ 



Uz 



k^A 

IT 



T 2z/5_ 



(62) 



In appendix |A] we show that the unique solution of ( !62l) that satisfies ( !6TI) is 
Uz = ikSut 



Ji(C)-fRe{exp(-iJC 



+ ^k'S'T, 



/2(C)-fRe^(7-if)exp(-^C 



(63) 



where 7 is the Euler-Mascheroni constant, 7 = 0.577..., ( = exp((zo — z)/6), as in §3.2.H 
and 



/i(C) 

/2(C) 



'sds 



e 

00 



sds 7 + In s 
1 + s^ 



(64) 
(65) 



Figure [TT] shows the vertical profiles of Uz and u^, and their analytical counterparts, from 
the same numerical solution shown in the top panel of figure [lOl 
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From fl63|) we deduce a relation between the constants Tt, Ut, and Wt, 



lC=o 



4 8 



(66) 



It can be verified that (!66|) holds, to reasonable accuracy, in each of the solutions presented 
in §4.11 Furthermore, in each case the term in (166|) involving Tt is found to be negligible, 
and so (l66ll reduces to the relation ( l49l) derived for the unstratified transition layer. The 
vertical mass flux within the tachocline is therefore tied to the differential rotation much as 
in the unstratified case described in ||3l 

We can also use f l63|) to quantify the role of the tachopause in the global-scale heat flow. 
Within the tachopause the thermal energy equation flTTl) becomes 



di2" 



(67) 



after making the boundary-layer approximation (l34l) . Therefore the change in the vertical 
temperature gradient across the tachopause is 



dT 

d^; 



2 = 20 — <5 



zo+5 



2o-<5 



(68) 



The transition between the weakly stratified and moderately stratified regimes occurs when 
( |68|) is of comparable magnitude to the temperature gradient within the tachocline (see 
figure [Hb). 



4-2.2. Global solution 

We can use the relations ( 166|) and (!68|) to construct an approximate analytical solution 
for the global flow. The details of the solution procedure are set out in appendix |Bl We find 
that the vertical velocity Wt in the tachocline is 

-i^icz 

= kd /1-;A 4 • l'^"' 



which reduces to equation fl50|) in the absence of stratification, rij-^ = 0. The dimensionless 
"geometrical" factors Gi, G2 and G3 are given by flB30l) - flB32l) in the weakly stratified 
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regime and by (lB33|l - flB35|) in the moderately stratified regime. We refer to tliese factors as 



geometrical because of their dependence on the tachocline thickness, D = h — zq. However, 
they also depend on /q, the thermal diffusivity enhancement factor defined in fl2Ul) . Moreover, 
G2 depends on the lengthscale d defined by ([SI]), and hence on the forcing timescale Tc- 

We can use (169|) to quantify the boundary between the unstratified and weakly stratified 
regimes identified in §4.11 This boundary is located where the term in the denominator of 
f l69p involving n^z becomes as large as the term involving 6. In general the geometrical factors 
are of order unity or smaller, and so the unstratified regime corresponds to 

^<k/S. (70) 

For the solutions in figures [7] and [8], which have E^, = 10~^, k = 2, and 6 = 0.02, this 
condition is equivalent to cr < 10~^, in good agreement with the numerical results. 

The boundary between the weakly and moderately stratified regimes occurs where the 
change in the temperature gradient across the tachopause, given by fl68l) . becomes compa- 
rable to the temperature gradient within the tachocline. As described in appendix |Bl this 
condition can be expressed approximately as 

1^ ^ yikS') . (71) 

For the solutions in figures [7] and [S] this condition is equivalent to a 0.025, which is 
consistent with (although somewhat smaller than) the value a ~ 0.1 found in the numerical 
solutions. 

The bottom panels of figures [7] and [9] verify that the analytical prediction ( 169|) for the 
vertical flow velocity in the tachocline agrees with the numerical results described earlier. 
Since viscosity was neglected in the derivation of (IMj) this good agreement demonstrates 
that the weakly stratified and moderately stratified regimes are both inviscid, that is, viscous 
forces do not play a significant role in the dynamics. In particular, the solutions do not have 
an Ekman layer at the radiative-convective interface , z = h. We believe that the existence 



of su c h an Ekman layer in some p revious studies (e.g. iGilman &: MieschI 120041 : iRiidiger et al. 



20051 : iGaraud fc Brummelll 120081 ) arises from the treatment of the interface in those studies. 

In all the results presented here, the change in buoyancy frequency n at the interface is 

smoothed over a lengthscale A, representing the depth of convective overshoot, that greatly 

1/2 

exceeds the Ekman length 6e = EJ . In the other studies just mentioned, the interface was 
modelled either as an upper boundary with a fixed differential rotation, or by a change in n 
over a lengthscale ^ 5e- Although the precise thickness of the Sun's overshoot layer is not 
known, it is surely thicker than the Ekman length (zz/fi©)^/^ ^ 30 m, and so we argue that 
there is no Ekman layer in the solar tachocline. 
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4.3. Physical interpretation and discussion of the stratified resuhs 

It is u seful to compare the expression for Wi given by fl69l) to the corresponding expression 



derived by lGBlOl for meridional circulations pumped between the outer and inner convective 
zones of young lithium-dip stars (see their equation (23)). As in their model, the magnitude of 
the circulation induced by the convective stresses is proportional to Ucz, a weighted average of 
the prescribed differential rotation in the convection zone. It is then moderated by whichever 
term in the denominator of (!69l) is largest. The first term is always of order unity or larger, 
and so Wt can never exceed Mcz in magnitude. The second term, involving 5, limits Wt to a 
flow rate which can be accommodated by the magnetic transition layer, as discussed in ^ 
this is an example of the "mechanical" constraint mentioned in 311 The third term, involving 
rirz, represents the thermal constraint described in ^ This constraint limits Wt to a flow 
rate for which the temperature perturbations created by the advection of the background 
stratification can be balanced by thermal diffusion, thereby maintaining thermal equilibrium. 

The expression for Wt given by fl59]) can be simplified considerably in the parameter 
regime appropriate to the solar interior. In the solar tachocline we have n^^/E^ ~ 8 x 10^^, 
and so the tachocline is not in the unstratified regime according to the criterion (!70|) unless 
the tachopause is unrealistically thin, 5 < 10^^'^. Therefore the denominator of ( |69|) is 
dominated by the thermal term. The weakly stratified and moderately stratified regimes, 
according to (|7T]) . correspond to 5 < 10^^ and 5 > 10^^ respectively. The solar tachocline 
could plausibly be in either of these regimes. Fortunately, it can be shown that Gi, G2 and 
(j3 follow similar scaling laws in either case, differing only by factors of order unity. We can 
therefore consider both regimes simultaneously. 

We expect the tachopause to be much thinner than the tachocline, 5 <^ which is 
itself much thinner than the solar radius, ^ 1. Under these conditions it can be shown 
that 

G1/G2 ~ UkDf (72) 
and Gi/Gg ~ D/6 > 1 (73) 

in both the weakly and moderately stratified regimes. The right-hand side of (172|l corre- 
sponds to the ratio of the thermal adjustment timescales in the tachocline and convection 
zone respectively. In the solar tachocline, thermal adjustment by radiative diffusion has 
a timescale of about 10^ years, whereas the thermal adjustment timescale in the convec- 
tion zone, estimated using mixing-length theory, is about 1 year. We therefore assume 
fo{kD)^ ^ 1, in which case G^ ^ G 2,G ?,. This same assumption was made implicitly in 



the models of ISpiegel fc Zahnl (119921 ) and iGough &: Mclntyrd (Il998l ). in both of which the 



top of the tachocline was assumed to be isothermal. Under this assumption we find that 
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Gi ~ (kD)^, and so fl69|) becomes 



2E, 



kD'^ 



wt ~ -iucz ^ , „o • (74) 



We note that Wt then depends only indirectly on the structure of the interior magnetic field, 
via the tachocline thickness D. The dependence of Wt on D can be understood physically 
as a consequence o f therm al equili brium and ther mal-wind balance, as first discussed by 



Gough &: Mclntyrd (Il998l ) (see also iMcIntyrd 120071 . p 194). Since the vertical shear across 



the tachocline is of order Ucz/D, thermal-wind balance implies that there must be a tem- 
perature perturbation T of order Ucz/{kD). Such a perturbation cannot persist within the 
convection zone, where heat is transported very efficiently by convective motions, but can 
persist within the tachocline, where heat transport is less efficient. In the tachocline, tem- 
perature perturbations diffuse at the rate Ei^/D'^, and thermal equilibrium therefore requires 
that n'^^'^t ~ ^K,T/D^ ~ {E,,/ D'^){ucz/{kD)), from which fITIl) follows immediately. Using 
order-of-magnitude estimates for the tachocline thickness, D ~ 0.01, and differential rota- 
tion, Mcz — 0.1, as well as n^z/E^ ~ 8 x 10^^, we find from (174|) that Wt ~ 10~^, dimensionally 
Wt ~ 10~^ cms~^. We note, however, that this result is rather sensitive to the tachocline 
thickness D, which is not well constrained by observations. 



Our model therefore recovers iGough fc Mclntyrd 's scaling for the meridional fiow within 



the tachocline, and goes further by clarifying the physical assumptions under which that 
scaling is derived, and by identifying the necessary conditions for those assumptions to hold. 
Under more general conditions, the meridional fiow strength Wt is given by fl^^ . and (1711) 
should be regarded instead as an approximate upper bound on Wt- The more general formula 
(169|) may be of relevance to the interiors of solar-type stars that are more weakly stratified, 
or that have thicker tachoclines. 

Of course, we must be careful when applying the results of our idealized model, with its 
artificially confined magnetic field, to the solar interior, or indeed to the interiors of other 
stars. However, since the formula for Wt given by ([71]) has no explicit dependence on the 
structure of the background magnetic field Bq, we argue that this result should hold for 
more general, less artificial field configurations than that considered here. Moreover, the 
physical processes that give rise to (I71|) will be present in any self-consistent model of the 
solar interior, and so (1711) ought to be a robust result under solar-like conditions. 

We now ask whether the meridional fiows predicted by (!74l) would be of sufficient mag- 
nitude to confine an interior magnetic field against outward diffusion. In our idealized model, 
the interior field Bq is confined by an artificially imposed downwelling Uq. A first approxi- 
mation to the nonlinear magnetic confinement problem can be obtained by setting Uq = \wt\. 
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Using ([7]), we can then relate the tachopause thickness S to the tachochne thickness D: 



_ 2Ek 



(75) 
(76) 



This expression is equivalent to equation (7) in iGough &: Mclntyrd ( 119981 ). 



We can now predict the strength of the magnetic field within the tachopause, Bt say , 
using equation f HSjl . Here, our model differs significantly from that of iGough fc Mclntyre . 
They assumed that thermal-wind balance would hold within the tachopause, as well as in 
the bulk of the tachocline, and as a result the thickness of the tachopause in their model, 
(5gm say, was tied to the strength of the stable stratification. They found that 



-"GM 



2E. 



1/6 



(77) 



where At = Aioc(-2o) = B^ / {AuripoVLo) is the Elsasser number within the tachopause. In our 
model, on the other hand, the Lorentz force within the tachopause overcomes thermal-wind 
balance, and we find that 

2 \i/2 



6 



k^At 



(78) 



(see §3.2. ip . Hence, whereas iGough fc Mclntyrd find a very strong dependence of B^ on 
D, namely Bt oc (see their equation (8)), we find by combining (175]) and ( 175]) that 
Bt oc D~^, and more precisely. 



B^ 8ul 



El 



47rpoi?|f^| nikW^ E^ 



(79) 



Taking /c = 2, n^^/E^ ~ 8 x lO^^^ and E^ ~ 3 x IQ-^^ we find 



B^ 



^ Bt 



^ 3 X 10-° — 
\0.1 



Ur 
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-3 

G. 



(80) 
(81) 



This estimate is rather higher than the ~ 1 G field predicted by iGough fc Mclntyrd . We note 
that torsional Alfvenic oscillations of a ~ 10 00 G field cou ld explain the 1.3 year oscillation 
detected in the tachocline's angular velocity (iGoughlboOoh . 
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4.4. Guidance for the selection of parameters in numerical models 

In recent years there have been several attempts to achieve magnetic field confinement in 
self-consistent, nonlinear, global numerical models of the solar int erior. Various approache s 



have been taken, including axisymmetric steady-state calculations (iGaraud fc Garaudll20081) 



as we l l as two-dimensional and three-dimensional time-de pendent simulations (IGaraud fc Rogers 



20071 : [Rogers fc MacGregorll201ll : IStrugarek et al.ll201ll ). However, none of these numerical 



models has so far obtained solutions that are satisfactorily close to solar observations, and 
none has achieved magnetic confinement over an extended range of la titudes. These failur es 



have led some to conclude that the magnetic confinement picture of iGough fc Mclntyrd is 
unworkable. However, the analytical results presented here suggest that magnetic confine- 
ment by meridional fiows can be achieved in the parameter regime appropriate to the solar 
tachocline. We can also use our results (a) to explain the lack of magnetic confinement in 
existing numerical models, and (b) to guide parameter selection for future models. 

As shown in §4.1| the amplitude of meridional fiows within the radiative interior de- 
creases monotonically as the stratification is increased. Moreover, in a time- dependent model, 
the timescale for the burrowing of the fiows, given by ([3]), increases with stratification. In 
the strongly stratified regime, with a > 1, the burrowing of meridional fiows is even slower 
than viscous diffusion. But if a < 1, as is the case in the tachocline, then the meridional 
fiows are approximately inviscid. When modelling the solar interior, it is common prac- 
tice to impose a rotation rate and stratification profile that are close to solar, in order to 
ensure a realistic separation between the dynamical timescales. The viscous, thermal, and 
magnetic diffusivities {u, k, rj) are then made as small as possible, subject to computational 
constraints. However, adopting this strategy does not guarantee that a < 1. Indeed, if the 
ratio of buoyancy and rotational frequencies n^^ is chosen to match the true tachocline value, 
then the condition cr < 1 requires the Prandtl number v/k to be very small, 

z//k<^~10"^ (82) 

^rz 

This condition is satisfied in the solar interior, but is beyond the reach of current numerical 
simulations, which instead typically have a Prandtl number much closer to unity, in order to 
reduce the numerical stiffness of the equations. This places such simulations in the strongly 
stratified regime in which (a) meridional fiow velocity decays exponentially with depth below 
the convection zone, and (b) viscous stresses make a significant dynamical contribution, even 
if the Ekman number is small, i.e. even if Ej, ^ 1. 

The numerical difiiculty is, in fact, rather easily avoided by imposing a weaker stratifi- 
cation, and thus allowing for a larger Prandtl number. For example, if = 10 then fl82l) 
requires only that u/k, < 10~^, which is readily achievable numerically. We then expect 



- 29 - 



the meridional flow velocity in the tachocline to scale as fTMj) . provided that the radiative- 
convective interface remains approximately isothermal (see discussion below fl72l) ). 



A further numerical difficultly is the need to resolve both the thin tachocline and the 
thinner tachopause. This difficultly can be alleviated by allowing the tachocline to be some- 
what thicker than is observed in the Sun@ According to (176!) . to obtain a tachopause of 
thickness 6 = 0.01 and a tachocline of thickness D = 0.1, for example, we require 

|^#^-10. (83) 

If we choose n^.^ = 10, as suggested above, as well as Ucz — 0.1 and k = 2, then this requires 
a diffusivity ratio tj/k = E^/E^ ~ 10~^. The strength of the interior magnetic field must 
then be chosen in accordance with ( [79|) . which requires 

^* - 50E. . (84) 



1/2 

Finally, we must ensure that the Ekman length E^, is smaller than S = 0.01. 

All of the constraints just described can be satisfied by choosing E^ = 10^'^, E^ = 10^^, 
and Ely = 10^^, as well as n^z = 10 and B^/^Anpo) = 0.05RqQq. These parameter values 
should be numerically a chievable; in f a ct, th e suggested values of E^, E^, and E^, are very 



similar to those used by lBrun &: ZahnI (j2006[ ). 



Discussion and Conclusion 



Following the failure of several recent attempts to recreate the iGough &: Mclntyre 



tachocline scenario in global numerical models, our r nain goal in this w ork was to iden- 



tify in what parameter regime, if any, the results of iGough &: Mclntyrd apply. For this 
purpose we created a model of the tachocline that is sufficiently simple to have analytical 
solutions yet, we believe, incorporates enough of the relevant dynamics to yield quantitative 
predictions. 

We have identified four distinct parameter regimes that occur in our results as the 
strength of the stratification is increased. For solar parameter values, the results lie in either 



■^The strength of the Sun's stratification increases rapidly with depth below the convection zone. For 
models that use a solar- like vertical profile of riiz, care must be taken to ensure that the condition a < 1 
holds throughout the tachocline. 
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the "weakly stratified" regime or the "moderately stratified" regime (see §4.1. 2p . In both 
of these regimes the strength of the meridional flow within the tachocline is determined by 
a combinati on of thermal equili brium and thermal-wind balance, and follows the scaling 



predicted by lGough fc Mclntyrd . With realistic tachocline parameters, the downwelling flow 
is of sufficient strength to confine an interior magnetic field across a thin boundary layer, 
and the thickness of the tachocline is related to the strength of the magnetic field by flHTI) . 

By contrast, all previous attempts to model the tachocline numerically have been per- 
formed in the "strongly" stratified regime, in which the burrowing of meridional fiows is 
significantly reduced by viscosity. We believe that this explains the lack of field confine- 
ment in those models. To remedy the problem, we suggest an alternative set of numerically 
achievable parameters, which we predict will yield results that are much more consistent 
with solar observations. 

This project was initiated during the ISIMA 2010 summer program, funded by the 
NSF CAREER grant 0847477, the France-Berkeley fund, the Institute of Geophysics and 
Planetary Physics, and the Center for the Origin, Dynamics and Evolution of Planets. We 
thank them for their support. P.G. and T.W. were also supported by NSF CAREER grant 
0847477, and J.M. was supported by NCAR and the Geophysical Turbulence Program. 



A. Analytical derivation of the transition-layer solution 

The magnetic transition layer is described by equations (I36l) - (l42|) . except that (!38|) is 
replaced by in cases with stratification. After eliminating p, ity, and by, these reduce to 

2^ = -pAb^e-'/' (Al) 
az 

= ifcK,e-^/^ + --^ + -^ A3 
az dz^ 



= M,e-'/^+l*i + :^. (A4) 

az az^ 



We now define a new variable C, = exp((2o — z)/6), with zq given by (I44p . and so the domain 
—oo < {z — Zq)/5 < +00 maps onto (X) > C > 0. When written in terms of equations 
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(IXT])-(IS1) become 

duz , I A 



dC 2C A;52 V 2 V ^ dC^ ^ dC 

2"^dC 



5 V 2 " dC^ 
After eliminating 6a; and fe^, we find 



dC=^ 

.d^M^ A;252Tt dw. 



ik5u, = -C-J7# (A9) 



and ik6^ = + (^ + , (MO) 

which we combine into a single equation for it^, 

This equation is equivalent to equation (162|) . but expressed in terms of (. We are interested 
in the solution of ( lAllI) that matches onto the flow in the weakly magnetic region above, 
and vanishes in the magnetically dominated region below. These matching conditions can 
be written as 



Ux ~ Ut + ^ikSTt In ( 



Uz -> Wt 

k ^ 



> as C^O (A12) 



and Uz ^ as ( oo. (A13) 

Conditions (lA12p correspond to , plus the condition that the magnetic field perturbation 
vanishes above the transition layer. 

After multiplying equation flAll|) by C, and integrating once, we find 

C'^. + = const. - ie6'T, InC . (AM) 
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Using flA9l) . we can write this as 



+ ikS ( - ) = const. - ^kH^Tt In C • (A15) 



For compatibility with (1A12|) the constant on the right-hand side must equal ikSut + hk'^S'^Tt. 



We can then use flA6P and (lASP to eliminate dux/d( and Uz from (lAlSp . which leads to 



1 /A d 



Ux-u,- \ik5T, InC = ■ (A16) 

The condition that — )■ as C ~^ implies that the right-hand side of f lA16p . and hence 
also the left-hand side, must be o{Q. That is, both sides of flA16p must vanish faster than C, 
as C — )■ 0. Similarly, the condition that — )■ as C — ^ implies that both sides of (lASP are 
o(l). So we can use f lAOp to express all four matching conditions in ( ]A12p as 



Equation flA14p . with the constant on the right-hand side now identified as ik6ut + 



^k'^6^Tt, becomes 



+ ^ = ik6u, 1 + Ik'd'T, 1^ . (A18) 
We can write the general solution as 



Uz = ik6u,h{0 + lk^6Xh{0 



+ ci exp (^C) + C2 exp (^c) + cs exp c) + c, exp (^c) (A19) 
where ci, C4 are arbitrary constants, and /i(C) and /2(C) are the integrals 

A(C)=r^^. (A20) 



e^' l + s 



00 



In equation (IA2ip 7 is the Euler-Mascheroni constant, 

POO 

7 = -/ e-^nsds = 0.577... (A22) 
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The values of the constants ci, C4 are fixed by the matching conditions flA13|l and 
(lA17p . To determine their values we need to consider the asymptotic behavior of Ji and l2- 
It can be shown that 



h 
h 

dC^ 

d^ 



and 



> as C 



(A23) 



Ii ~ 
I1 ~ 



as C 00 . 



(A24) 



Matching condition flA13l) . together with flA24p . implies that Ci = C3 = 0. The remaining 
matching conditions flA17p . together with flA23p . then imply that 



Wt = if fc^Mt + f 7^ ^ Tt + C2 + C4, 
= i^k6u, + ^(7 + f )A;2(52Tt + cse^/^ + C4e-'-/^ 



and 



= i^fc^Mt + -5=(7 - j)k''6'% + C2e-'-/^ + c^e 



^i7r/4 



(A25) 
(A26) 
(A27) 



2V2 * 4^2^ 

Solving the three equations ( 1A25|) - (]A27I) fixes the values of C2 and C4, and also imposes a 
condition on Wt, 

Wt = -if Mut - f 7fc'5'Tt . (A28) 

The solution for -u, is 



ik6ut 



/i(C)-fRe^xp(-i^C 



/2(C)-fReU7-if)exp(-iJC 



(A29) 



We can also calculate the change in the temperature gradient across the transition layer. 
This is given approximately by fl55]) . and more precisely by 



dT 
dJ 



z=+co 



z=+oo 



Wt dz 



dC 



Uz—r - Wt — 



((f)^iA:Vr,-7^t) . 



(A30) 

(A31) 
(A32) 
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B. The global solution, and the vertical flow velocity in the tachocline 



As in the work of iGBlOl . we construct an approximate global solution of fl22|) - fl29|) by 
finding the general solution in each region of the domain and then matching these solutions 
across the boundaries. In our case the boundaries between the regions are aX z = h and 
z = zq^ and are known a priori, with h given by the background stratification and Zq given 
by dm). 



B.l. The general solutions in each region 
B.1.1. Solution in the convection zone, z G [/i, 1] 
In the convection zone the governing equations are well approximated by 



-2Uy = 




(Bl) 




2ux = 


—ikp 

rc 


(B2) 


= 


dp - Uz 
dz Tc 


(B3) 


= 


ikUy + — — 
dz 


(B4) 


= 




(B5) 



The general solution for the temperature perturbation can be written as 

T = a cosh k{z — 1) + b sinh k{z — 1) 



(B6) 



where a and b are integration constants. From the boundary condition T = at 2; = 1 we 
deduce immediately that a = 0. Combining the remaining equations yields 



[d^ d^^""' 



-T + 



2i d-Ucz 
k dz 



rc_ 
d? 



(B7) 



where d is the lengthscale defined in equation (|5T1) . We write the general solution for u-^ as 



2i7, 

' kd^ 



dz' Ucz{z') cosh 



z — z 

d 



sinh k{z—l)+A cosh 

4r, 



z-l 
d 



+B sinh 



z-l 



d 
(B8) 

where A and B are two additional integration constants. The boundary condition w = at 
z = 1 implies that A = 0. (Since we have neglected the viscous terms in ( ]Bl[) -( lB3l) we cannot 
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impose the stress- free boundary condition at 2; = 1. Including the viscous terms would lead 
to an Ekman-type boundary layer forming at z = 1, but the effect on the solution within 
the bulk of the convection zone would be of order <^ 1.) 

Finally, the pressure perturbation is found to be 

p=^'-l<^ (B9) 

2i f\ , , , fz'-z\ hkd? , , Bd f z - l\ , 

= — — / dz Ucz\z ) smh — ; — + cosh kiz — 1) cosh — ; — . (BIO) 

kdj^ \ d J At^ Tc \ d J 



B.1.2. Solution in the tachocline, z G [zQ,h] 
Within the tachocline we have 

-2uy = (Bll) 
2u^ = -ikp (B12) 

= -^ + f (B13) 

az 



= ikuy + ^ (B14) 
dz 



d^T 



(IB11|) and (]B14p together imply that Uz is a constant, Uz = Wt say. The remaining equations 



can then be integrated, and the general solution written as 

Uz = Wt 



T = \^, + ^w,jcos\ik{z-Zo) + KsiYi\ik{z-Zo)-^Wt 
Ux = Mt - li ( Tt + ji^Wt sinh k{z - 2:0) - iiJ^(cosh k{z - zq) - 1) + ^'l, ^ {z - zq) 



^ = + i + ^^y;^^ sinh k{z - Zq) + y (cosh k{z - zq) - 1) - ^^(^ " ^0) 

where Ut and Tt are the values of Ux and T at the bottom of the tachocline, and K is an 
additional integration constant. 
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B.1.3. Solution in the magnetically- dominated region, z < zq 

In this region the differential rotation and meridional ffow both vanish. The temperature 
perturbation T therefore satisfies 

d^f 

= ^-eT (B16) 
T = a cosh kz + P sinh kz . (B17) 

Since T = at the lower boundary z = 0, we must have a = 0. 



B.2. Matching conditions 

The values of the seven integration constants b, B, wt, Ut, Tt, K and /3 are now de- 
termined by applying matching conditions across the interfaces z = and z = h. At the 
radiative-convective interface, z = h, we impose that Uz, p, T, and f{z)^ are all contin- 
uouslll These are the same continuity conditions imposed by GBIO . except that we use a 



more realistic thermal energy equation ([9]), and as a result we require the continuity of the 
heat flux, rather than the temperature gradient. Across the magnetic transition layer, at 
z = Zq, we impose continuity of T and the relations (IA28P and f lA32p derived from our 
transition-layer solution. 

The matching conditions lead to the following seven relations between the integration 



* By imposing that iiz is continuous a.t z — h we neglect any gyroscopic pumping within the overshoot 
region. This is reasonable provided that the overshoot depth A is not too large. 
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constants: 
2ir, 



kd? 



dz' Ucz{z') cosh 
= Wt 



+ - — sinh k{l — h) — B sinh 



1 - h 
d 



d 



bkd^ , Bd , 

H — — r cosh k[l — h) cosh 

4r2 



1 - h 
d 



2iMt 1 
— - + - 
k k 



nt 



K 



n 



6 sinh k{l -h)= \ Tt + y^^t cosh kD + Ksinh kD 



Tt + T^u't sinh + —(cosh kD - I) - 
k-^E^ 



k^E., 



bk cosh A;(l — h) 
Tt 



k 
Jo 



) sinh kD + —K cosh kD 
Jo 



(3 sinh kzo 

Kk = Pkcoshkzo + ^ ((f )' l^'^^Tt - 7u;t) 



.[EkSut - llk'5'Tt 



(B18) 

D 
(B19) 
(B20) 

(B21) 
(B22) 
(B23) 
(B24) 



where D = h — zq is the tachochne thickness. 



We now seek an exphcit expression for the value of the vertical flow Wt within the 
tachocline. We begin by eliminating B between (IBlSp and (IBlOp . which leads to 



ikd 

— Wt coth 



1-h 
d 



+ 



in^w^ 



FE. 



nz 



kD-i[Tt + jy^wt sinh kD - iir(cosh kD - 1) 



2ucz - 2ut + 



ibkd 
4^ 



sinh k{l — h) 
tanh(i^) 



FE« 



— fcficosh — h) 



(B25) 



where Ucz is the weighted average of the forcing in the convection zone defined by equation 
(l52l) . The remaining six integration constants, including Wt, therefore depend on Ucz{z) only 
through its average Ucz- Hence the entire solution below the convection zone is determined 
by Ucz- 

Next, we combine flB22p and flB23p to eliminate /3. The result can then be combined 
with ( lB20p and ( ]B2ip to express 6, Tt and K in terms of Wt- The general result is rather 
complicated, so for simplicity we describe only two limiting cases, which correspond to the 
"weakly stratified" and "moderately stratified" regimes identified in §4.1.2[ 



- 38 - 



B.2.1. The weakly stratified regime, ra^^/E/t ^ l/{k6^ 



In this regime we find that the term in (IB23P involving n-rz becomes neghgible. Since 



this term arises from the change in the temperature gradient across the transition layer, this 
regime corresponds to the "weakly stratified" regime described in §4.1.2[ After neglecting 
this term, we find that 



(B26) 



^ n^^Wt 1 — sech k{zo + D) cosh kz^ 

/c^Ek /o cosh — h) tanh k{zQ + D) + sinh k{l — h) 
^ n^^Wt /o tanh kzQ[l — sech kD] + tanh — h) tanh kz^ tanh kD (B27) 

/c^Ek /o [tanh kzQ + tanh kD] + tanh k{l — h) [tanh kzQ tanh kD + 1] 
^ n^^Wt /o[l — sech kD] + tanh k{l — h) tanh kD (B28) 

/c^Ek /o [tanh /c^Tq + tanh kD] + tanh — h) [tanh /c2;o tanh /cD + 1] 



B.2.2. The moderately stratified regime, n'^^/E^ ^ l/{k6^) 



In this regime the transition-layer temperature Tt turns out to be smaller than the 
value given by (IB27P by a factor E^/{n'^^k6^) <^ 1. To good approximation, therefore, we 
can neglect Tt in each of the matching conditions ( IB18p - flB24p . which is equivalent to setting 
= in (!B26l) -( !B28l) . 



After eliminating ut between equations ( 1B25P and ( 1B24I) , and applying the formulae for 
b, Tt and K just derived, we arrive arrive at an explicit formula for Wt, 



where 



Go 



-lUc 



Wt 



kd , 
— coth 

2Tr 



h 



d 



nz 



7ik5 2A;2E, 



- [Gi + G2 ~ G3 



Gi = kD- 



/o[2 — 2sech kD + tanh kzQ tanh kD] + tanh k{l — h) tanh kD 
/o[tanh kzQ + tanh kD] + tanh k{l — h) [tanh kzo tanh kD + 1] 



kd 



k^d^ 



kd — 



tanh ^(1 — h) 
tanh(i-^) 



sech kh cosh kzo 



/o tanh kh + tanh k(l — h) 



G3 = ik5 



/o[l — sech kD] tanh kzQ + tanh A; (1 — h) tanh kD tanh kzQ 
/o[tanh kzQ + tanh kD] + tanh k{l — h) [tanh kzo tanh kD + 1] 



(B29) 



(B30) 
(B31) 
(B32) 
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in the weakly stratified regime, and 

/o[2 - 2sech kD] + tanh k{l - h) tanh kD 



Gi = kD 



fo tanh kD + tanh k{l — h) 



kd 



kd 



k^d^- 
G3-O 

in the moderately stratified regime. 



tanh k{l — h) 
tanh(i^) 



sech kD 



fo tanh kD + tanh k{l — h) 



(B33) 

(B34) 
(B35) 
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Fig. 4. — Profiles of vertical velocity for a series of unstratified solutions with varying magnetic field 
strengths. The Elsasser number A was increased from 10^ to 10"^^ in multiplicative increments of 10^. Other 
parameters are as in figure [S] In each solution, the magnetic transition layer is located where the local 
Elsasser number Aioc = Ae'^^/"^ ~ 2/{k6)'^ (see pXTj) . 
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Fig. 5. — The vertical profiles of iiz and Ux within the magnetic transition layer. The solid lines correspond 
to the numerical solution shown in figure|3l the symbols correspond to the analytical boundary-layer solution. 
The Uz profile has been multiplied by 10 to make the two profiles visible on the same scale. The vertical 
axis represents the boundary-layer coordinate [z — zq) and so the convection zone lies far above the region 
plotted. 
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Fig. 6. — Top panel: Vertical velocity profiles for various values of the magnetic scaleheight 5. Other 
parameters are the same as in figure [3l The dotted lines show the asymptotic behavior (|48|) of the boundary- 
layer solution (|46| in the magnetically dominated region. Bottom panel: The vertical velocity Wt in the 
weakly magnetic region predicted by (|50p . The symbols show the value of at z — 0.65 in each of the 
numerical solutions shown in the top panel. 
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Fig. 7. — Top panel: Profiles of Uz{z) from a series of solutions with varying a. All the solutions have 
/o ~ 10^ in (pn|) . and the other parameters are the same as in figure |31 We use solid lines for the unstratified 
solutions (cr = 10^^, lO^"*, lO""^), dashed lines for the weakly and moderately stratified solutions (cr = 
10-2■^10-^10-^■^10-^10"°■^), and dotted lines for the strongly stratified solutions (cr = 10°,10\l02). 
Bottom panel: The dotted and dashed lines show the vertical velocity Wt in the tachocline predicted by 
for the weakly stratified and moderately stratified regimes respectively; a solid line is used to indicate the 
range of a values in which each is expected to apply. The symbols show the value of Uz &i z = 0.5. The 
vertical lines indicate the boundaries between the unstratified, weakly stratified, moderately stratified and 
strongly stratified regimes. 
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Fig. 8. — Profiles of Ux{z) and T{z) from the same solutions shown in flgure[7l The linestyles are the same 
as in figure [T] Note that the vertical axes are linear. 
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Fig. 9. — Top panel: Profiles of Uz{z) from a series of simulations with A increasing from A = 10^ to 
10^* in multiplicative increments of 10^, for a = 10~^'^ (left column) and a = lO"*^'^ (right column). Other 
parameters are as in figure [T] Bottom panel: The vertical velocity Wt in the tachocline predicted by (j69p 
in the weakly stratified regime (left column) and moderately stratified regime (right column). The symbols 
show the value of at z = 0.65, just below the base of the convection zone, for each of the profiles in the 
top panel. 
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Fig. 10. — Meridional cross-sections through the solutions in figure [7] with cr = 10~^'^ and cr = 10^, 
corresponding to the weakly and strongly stratified regimes respectively. The left column shows the azimuthal 
velocity and meridional streamlines. The right column shows the azimuthal magnetic field and poloidal field 
lines. Different streamlines have been plotted in the two cases because in the strongly stratified regime the 
meridional fiow within the tachocline is weaker by many orders of magnitude. 
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Fig. 11. — The vertical profiles of iiz and Ux within the tachopause. The solid lines correspond to the 
numerical solution shown in the top panel of figure I10( the symbols correspond to the analytical boundary- 
layer solution. The Uz profile has been multiplied by 50 to make the two profiles visible on the same scale. 



